1 Disclaimer

The following work is preliminary and subject to ongoing investigation.

These analyses do not necessarily reflect the point of view of ICCAT or the FAO and in no way anticipate ICCAT or FAO future policy.

2 Objective

Develop a flexible and powerful indicator of stock status for tunas, billfish and sharks that can be tested theoretically and empirically.

3 Project details

Title: ‘Simulation Testing Ecosystem Indicators: Support to ICCAT’s Ecosystem Approach to Fisheries Management (EAFM)’

Term Dec 2023 - June 2026
Funding body International Commission for the Conservation of Atlantic Tunas (ICCAT)
Funding stream FAO Global Environmental Facility (GEF) Common Oceans ABNJ Tuna Project Fund
Call No. # 11903 / 2023
Project Partners Blue Matter Science Ltd.
Blue Matter Team Tom Carruthers, Adrian Hordyk, Quang Huynh
ICCAT Collaborators Nathan Taylor

4 Introduction

To meet the requirements of the precautionary approach and the ecosystem approach to fisheries management (EBFM), indicators of stock status are needed for secondary species, defined here as those lack sufficient data or capacity to conduct routine stock assessments (see Carruthers et al. 2024 for a detailed problem statement). These indicators must be theoretically sound and be validated empirically.

The EcoTest project aims to use simulation modelling to identify data and algorithms that can inform stock status of secondary species and then validate these indicators empirically in cases where there have been defensible stock assessments.

A well-documented, defensible and transparent framework is needed to support tactical decision making to move beyond the single species assessment paradigm and make progress towards the essential goals of EBFM.

Previous work synthesized the dynamics of six stock assessments and consolidated those into a multi-species, multi-fleet operating model (Huynh et al. 2022). Those operating models were used to simulate various future scenarios for fishing and stock dynamics to develop Indicator 1: an indicator system specific to the longline case study that could predict stock status for the 6 species without undergoing a full stock assessment (Carruthers et al 2024).

The key advance of this earlier work was to use deep learning to train neural networks on the posterior-predicted data of a large number of simulations (a training dataset). The ETest package is built around an evoluation of that approach, using simulated posterior predicted data from a much more general set of simulations that include population and fishing dynamics for a wide range of pelagic highly migratory species that are caught in Atlantic fisheries.

These indicators focus on the prediction of stock status relative to a productive level expressed as spawning stock biomass relative to spawning stock biomass corresponding to MSY levels. This is a key reference level for evaluating stock status for the majority of stock assessed by ICCAT. Expert working groups also confirmed that there is interest in developing indicators that can accommodate various numbers of fleets and stocks, and that these should also apply to varying data availability scenarios with respect to length, age and relative abundance data.

5 Overview of the EcoTest approach

5.1 Propose a suite of data inputs

EcoTest can accept a range of data inputs (features) that may be available for secondary species derived from data streams such as total annual catches, CPUE indices, mean length in catches and mean age in catches. These can be provide along with life-history (e.g. somatic growth, natural survival) and fleet-specific (e.g. selectivity) characteristics to potentially improve reliability in interpreting these data inputs.

Where multiple fleets can be identified, data inputs can be provided for up to 3 individual fleets.

Where multiple stocks can be assumed to have been exploited by a common fishery with comparable exploitation history, multi-stock data inputs can be used including catch ratios of one stock to another.

5.2 Format data inputs for use with EcoTest neural networks

A training dataset was generated from a large number of multi-stock, multi-fleet simulations spanning a range of life-histories and fishery dynamics typical for tuna, billfish and sharks. In each simulation more than 400 data inputs were calculated along with the known, true simulated stock status (spawning biomass relative to spawning biomass at MSY levels). The data inputs relate to trajectories, levels and variability of data such as catches, CPUE indices and mean length in catches and are processed and formatted in such a way that they can be calculated for any fishery (e.g. mean lengthin catches is phrased as a percentage of asumptotic length, slopes in data are slopes in log data etc).

Some data inputs relate to data from fishing fleets, some relate to fishery characteristics, and others relate to the life hsitory of the stock(s) (see Table 1).

5.3 Train and artificial neural network for the specific suite of data inputs selected by the user

The large training dataset is subsetted to include only those data inputs provided by the user. An artificial neural network is then trained on the substted training dataset.

During training a validation dataset is used to check for model overparameterization and confirm predictive ability for independent data not used by the training algorithm.

After training, the expected indicator performance is evaluated according to a completely independent testing dataset.

This step can be seen as (a) a theoretical check that there is information in the features (input data) provided to estimate stock status and (b) a quantitative theoretical evaluation of the expected performance of the indicator.

A result of this step is a fitted model and also a standardized version of the user’s dataset that has the same standardization process as the training dataset.

5.4 Estimate stock status using the indicator

The product of step 3 (fitted model and standardized version of the dataset) can be used to estimate stock status.

If a historical trend in empirically estimated stock status is provided, the estimates of the indicator can be compared with other empirical estimates such as those arising from stock assessments.

6 Installation

6.1 R and RStudio

The Keras, Tensorflow and Miniconda libraries on which EcoTest methods are based are highly sensitive to the version of R and RStudio.

It is therefore highly recommended that users install the very latest version of R and RStudio before continuing (you will almost certainly have difficulties otherwise).

R can be downloaded from the R Project for Statistical Computing webpage.

Currently the ETest package is compatible with R version R version 4.5.1 (2025-06-13 ucrt).

RStudio can be downloaded from the Posit webpage.

Currently the ETest package is working in RStudio version 2025.05.1 Build 513.

The ETest package has only been tested on a Windows x64 system.

6.2 Installing R Packages

Open RStudio and run the following lines of code at the R command prompt:

# Keras, Tensorflow, Miniconda
install.packages("keras3")
library(keras3)
install_keras()

# EcoTest package ETest
install.packages("remotes")
remotes::install_github('blue-matter/ETest')

6.3 Test Installation

To test your installation run the following code

library(ETest)
test = train_ind(Shark_1_data)

If your packages are working you should see something like this:

7 Worked Examples

7.1 Quick Start

Lets just assume we have some data on catches, recent CPUE and a bit of life history info. Remember, these names have to match the EcoTest codes for indicator features (see Table 1 below).

library(ETest)

dat = data.frame(
  K_s1 = 0.25,         # von B. somatic growth (stock 1)
  M_K_s1 = 1.05,       # Instantaneous natural Mort. / von B. K (stock 1)
  L50_Linf_s1 = 0.74,  # Length at 50% maturity (spawning fraction) / asymptotic length (L-infinity) (stock 1)
  maxa_s1 = 11,        # Age at 5% survival (stock 1)
  L5_L50_s1_f1 = 0.2,  # Shortest length at 5% selectivity relative to length at 50% maturity (stock 1)
  C_rel_s1_f1 = 0.67,  # Current catches relative to historical average (stock 1, fleet 1)
  C_g10_s1_f1 = 0.01,  # Slope in log catches over the last 10 years (stock 1, fleet 1)
  I_rel_s1_f1 = 0.43,  # Current CPUE relative to historical average (stock 1, fleet 1)
  I_g20_s1_f1 = -0.09, # Slope in log CPUE over the last 20 years (stock 1, fleet 1)
  ML_rel_s1_f1 = 0.9   # Current mean length in catches relative to historical average (stock 1, fleet 1)
)

input = list(dat=dat)
my_indicator = train_ind(input)

You can see that the performance of this indicator isn’t stellar however it correctly classifies stocks below half SSBMSY (an ICCAT limit reference point) and above SSBMSY (an ICCAT target reference point) in about 60% of simulations.

For now we are going to assume we find this acceptable (you might well not!).

We can now obtain a point estimate of stock status using our indicator and the dataset provided:

my_pred = pred_ind(my_indicator)
my_pred

Since only one set of data were provided there is only a point estimate of stock status. However most datasets can be sampled to provide stochastic inputs (bootstrapping etc).

For the purposes of demonstration lets just great some samples of data so you can plot indicator estimates that include uncertainty:

nsamp = 50                                                  # 50 stochastic draws of the dataset
nobs = ncol(dat)                                            # data set has nobs features (data types)
log_norm_err = rlnorm(nobs*nsamp,0,0.2)                     # demo error 
data_repped = rep(as.numeric(dat),each=nsamp)               # duplicate data for nsamp rows
dat_s = array(data_repped * log_norm_err,c(nsamp,nobs))     # combine data with errors
dat_s = as.data.frame(dat_s)                                # make into data frame
names(dat_s) = names(dat)                                   # make sure to copy over labels

input_s = list(dat=dat_s)                                   # make an input object where position one in the list is the data
my_indicator_s = train_ind(input_s)                         # train the indicator
my_pred_s = pred_ind(my_indicator_s)                        # make a status prediction for each of the 50 samples of data
hist(my_pred_s, xlab="Estimated SSB/SSBMSY", main="Demonstration Stochastic Indicator", col="cornflowerblue", border="white")

 

7.2 Get a point estimate of stock status from single species data

In this example we will process some example data, train an indicator and estimate stock status using it.

The ETest package comes with an example dataset ‘Some_data’ for these worked examples. The dataset has catches, CPUE indices, mean length, life history and selectivity information for three fleets simultaneously fishing three fleets. We are not necessarily going to assume we have all these data, but they are there to experiment with.

To make creating time series data features easier, ETest has a function ts_features() which extracts the levels and slopes processed in the same way as the training dataset.

stock = 1                                # Lets just try this for one stock and ...
fleet = 1                                #    ... one fleet. 
Catch = Some_data$Catch[stock,fleet,]    # Catch data is a 3D array [stock, fleet, year]
plot(as.numeric(names(Catch)),Catch)     # what do these catch data look like?
Catch_f = ts_features(Catch, "C")        # lets capture the rel, g5, g10, g20 and g40 features
Catch_f

Index = Some_data$Index[stock,fleet,]    # CPUE index data is a 3D array [stock, fleet, year]
Index_f = ts_features(Index, "I")        # captures the rel, g5, g10, g20 and g40 features

ML = Some_data$Mean_length[stock,fleet,] # Mean length data
ML_f = ts_features(ML,"ML")              # captures the rel, g5, g10, g20 and g40 features

fleet_data = as.data.frame(c(Catch_f, Index_f, ML_f))                  # combine fleet data
names(fleet_data) = paste0(names(fleet_data),"_s",stock,"_f",fleet)    # label correctly 

LH_sel_data = data.frame(K_s1 = Some_data$K[stock],                    # combine life history and selectivity info
                       M_K_s1 = Some_data$M[stock]/Some_data$K[stock], # ratio of M to K 
                       maxa_s1 = Some_data$maxa[stock],                # age at 5% natural survival
                       L5_L50_s1_f1 = Some_data$L5[stock,fleet]/Some_data$L50[stock],   # shortest length at 5% selectivity relative to L50
                       LFS_L50_s1_f1 = Some_data$LFS[stock,fleet]/Some_data$L50[stock]) # length at 95% selectivity relative to L50

input2 = list(data=cbind(LH_sel_data, fleet_data))                     # Combine inputs into a single input object

Now that we have our dataset, extracted from three time series (Catch, a CPUE index and Mean Length) in addition to selectivity and life history parameters, we can set up the indicator:

my_indicator2 = train_ind(input2)  # train indicator, standardize input data

And then get an estimate of stock status:

pred_ind(my_indicator2)  # estimate stock status from trained indicator and submitted dataset

7.3 Conduct empirical evaluation using ETest datasets

The ETest package comes with a number of real fishery datasets loaded. You can list these using the R function data():

data(package="ETest")

There is an object TD which is the full training dataset. You can always look at that if you want to check the magnitude and range of data inputs used for training the neural networks.

The other objects are labelled _io, _data and _retro.

_io objects such as Shark_1_io are list objects that contain all of the inputs and outputs to a Stock Synthesis 3 assessment.

_data objects are list objects that contain ETest formatted features obtained from those _io stock synthesis runs (position 1, ‘inputs’) and also the stock synthesis estimated stock status (position 2, ‘Brel’).

_retro objects are a list of _data objects npeels long which strip back the data from the stock synthesis runs to allow retrospective behavior of indicators to be evaluated.

It is simple to conduct indicator training and prediction, including retrospective analysis:

# Single stock status estimate using data to most recent year:
Shark_1_ind = train_ind(Shark_1_data)
Shark_1_pred = pred_ind(Shark_1_ind)
hist(Shark_1_pred)

# Retropsective stock status estimate using npeels of the SS3 data:
Shark_1_ind_retro = train_ind(Shark_1_retro)
Shark_1_pred_retro = pred_ind(Shark_1_ind_retro)

7.4 Use an existing stock synthesis assessment to conduct an empirical test

To obtain all relevant outputs for a current Stock Synthesis model you will a directory where an assessment has been run and you will need to install the latest version of the r4ss package.

wm_io = all_ss3("C:/myassessments/ATL_white_marlin_2025")

Once you have the assessment ripped into R, you need to pick what fleets and indices you wish to extract. You can use the ss_names() function to list the fleets / surveys. You then have to select the relevant fleets and corresponding cpue indices and then make your formatted EcoTest indicator features using SS_2_ET():

ss_names(wm_io)
wm_data = SS_2_ET(wm_io, Fnam = c("JPN_LL_N","US_PL_N"), Inam = c("Surv_JPN_LL_N","Surv_US_PL_N"))

The wm_data object now has features for three fleet data sets: the Japanese longline north the US pole and line north and all other fleets combined. Selectivities for all other fleets combined are weighted by historical average catches.

The function SS_2_ET() samples from time series and parameters to provide a stochastic feature dataset in order to obtain variance in stock status estimates.

You can now run the indicator for that dataset:

wm_ind = train_ind(wm_data)
wm_pred = pred_ind(wm_ind)

or alternatively you can do the same steps but conduct a retrospective analysis using the function SS_2_ET_Retro():

wm_retro = SS_2_ET_Retro(wm_io, Fnam = c("JPN_LL_N","US_PL_N"), Inam = c("Surv_JPN_LL_N","Surv_US_PL_N"),npeels=8)
wm_ind_retro = train_ind(wm_retro)
wm_pred_retro = pred_ind(wm_ind_retro)

7.5 Multi-stock indicators

In general, fisheries stock assessment has failed to incorporate information across multiple species in the estimation of stock status. This is a missed opportunity.

For example if the denominator of CPUE indcies of abundance (effort) is comparable in trend among two species, their catch ratio informs their relative depletion (e.g. if the catch ratio of stock 1 : stock 2 has declined by 50%, then this infers that stock 1 is at half the vulnerable stock size as stock 2 (assuming you standardized your data correctly).

This means that depletion on one stock may be informed by data informing depletion on the other stock. As the number of stocks increases (that adhere to this assumption of approximately comparable effort trends) an increasibly narrow set of depletion conditions may prevail. To demonstrate this lets see how well only catch data can inform stock status across three species:

fleet = 1                                #    ... one fleet. 
Catch_f1 = ts_features(Some_data$Catch[1,fleet,],"C")
Catch_f2 = ts_features(Some_data$Catch[2,fleet,],"C")
Catch_f3 = ts_features(Some_data$Catch[3,fleet,],"C")
cvec = c(unlist(Catch_f1),unlist(Catch_f2),unlist(Catch_f3))
msdat = data.frame(matrix(cvec,1))
names(msdat) = paste0(names(cvec),"_s",rep(1:3,each=5),"_f1")
input = list(dat=msdat)
catch_ind = train_ind(input)

Since there is nothing providing an indication of absolute depletion the indicator based on catches alone performs very poorly.

fleet = 1   #    ... one fleet. 


Index_s1 = ts_features(Some_data$Index[1,fleet,],"I")
Index_s2 = ts_features(Some_data$Index[2,fleet,],"I")
Index_s3 = ts_features(Some_data$Index[3,fleet,],"I")
Index_rels = data.frame(I_rel_s1_f1 = Index_s1$I_rel,
                        I_rel_s2_f1 = Index_s2$I_rel,
                        I_rel_s3_f1 = Index_s3$I_rel)
wIdat = list(dat = cbind(msdat, Index_rels))
catch_wI_ind = train_ind(wIdat)

Adding just the current index levels relative to historical means, tidied the indicator up somewhat.

Now lets find out how much of that was just the index levels:

catch_onlyI_ind = train_ind(list(dat=Index_rels))

As expected, the relative index levels were doing the majority (but not all) of the lifting.

8 Summary of EcoTest Indicator Features (data types)

Table 1. Indicator Features (data types)

Feature Type Description
K Stock von Bertalanffy growth parameter (annual)
M_K Stock The ratio of natural mortality rate to von B. K
maxa Stock Maximum age. The age corresponding to 5% natural survival (the age before which 95% of individuals have died from natural mortality)
L50_Linf Stock The ratio of length at 50% maturity (spawning fraction) to asymptotic length (von B. L-infinity)
L5 Fleet The shortest length at 5% selectivity phrased as a % of asymptotic length (L-infinity)
LFS Fleet The shortest length at 95% selectivity phrased as a % of asymptotic length (L-infinity)
VML Fleet The vulnerability at maximum length (‘dome-shapedness’)
I_rel Data The ratio of current index level to mean historical index level
I_g5 Data The slope in log index over the most recent 5 years
I_g10 Data The slope in log index over the most recent 10 years
I_g20 Data The slope in log index over the most recent 20 years
Isd1 Data Mean absolute residual error in a log CPUE index from a loess smoother with effective number of parameters of 40% index length
Isd2 Data As Isd1 but ENP of 20%
Isd3 Data As Isd1 but ENP of 10%
C_rel Data The ratio of current catches relative to mean historical levels
C_g5 Data The slope in log catches over the most recent 5 years
C_g10 Data The slope in log catches over the most recent 10 years
C_g20 Data The slope in log catches over the most recent 20 years
Csd1 Data Standard deviation in residual error around log catches from a loess smoother with ENP 20% catch series length
CF Data The fraction of total historical catches (taken by this fleet).
ML_rel Data The ratio of current mean length in catches relative to historical average mean length in catches. Mean length is phrased as a % of asymptotic length (L-infinity)
ML_g5 Data The slope in log of mean length of catches over the most recent 5 years
ML_g10 Data The slope in log of mean length of catches over the most recent 10 years
ML_g20 Data The slope in log mean length of catches over the most recent 20 years
ML_Linf Data The ratio of current mean length in catches relative to asymptotic length (von B. L-infinity)
ML_L50 Data The ratio of current mean length in cathes to length at 50% maturity (spawning fraction)
MA_rel Data The ratio of current mean age in catch relative to mean historical levels
MA_g5 Data The slope in log mean age in catch over the most recent 5 years
MA_g10 Data The slope in log mean age in catch over the most recent 10 years
MA_g20 Data The slope in log mean age in catch over the most recent 20 years
MV_rel Data The standard deviation in observed caught lengths in the current year relative to the historical average.
MV_g5 Data The slope in log standard deviation in catches over the most recent 5 years
MV_g10 Data The slope in log standard deviation in catches over the most recent 10 years
MV_g20 Data The slope in log standard deviation in catches over the most recent 20 years
FM_rel Data The ratio of current fraction mature fish in catches relative to historical mean fraction of mature fish in catches
MV_g5 Data The slope in log of fraction mature fish in catches over the most recent 5 years
MV_g10 Data The slope in log of fraction mature fish in catches over the most recent 10 years
MV_g20 Data The slope in log fraction mature fish in catches over the most recent 20 years
CR Data The current catch ratio (first stock divided by second stock). R(y) = C(1,y) / C(2,y). CR is R from the most recent year, R(ny)
CR_mu Data The catch ratio (first stock divided by second stock) over the first 20 years R[1:20]
CR_rel Data CR / CR_mu
CR_s5 Data The slope in log CR over the most recent 5 years
CR_s10 Data The slope in log CR over the most recent 10 years
CC20 Data Correlation of the in residual fit of a loess smoother to log catches with ENP = 20% among stocks over the most recent 20 years.
CC40 Data As CC20 but over historical years between 21 and 40 years ago.
CC60 Data As CC20 but over historical years between 41 and 60 years ago.

9 Resources

The openMSE website contains extensive documentation about openMSE and approaches such as RCM.

13 References

Carruthers, T.R. 2018. A multispecies catch-ratio estimator of relative stock depletion. Fisheries Research. 197: 25-33, https://doi.org/10.1016/j.fishres.2017.09.017.

Carruthers, T.R., & Taylor, N.G. 2025. EcoTest Indicator 2: A general purpose stock status indicator for sharks, billfish and tunas.

Carruthers, T.R., Huynh, Q., Taylor, N.G. 2024. EcoTest Phase III: Simulation testing Ecosystem Indicators. SCRS/2023/087.

Huynh, Q., 2025. Rapid Conditioning model. Available from https://openmse.com/tutorial-rcm/

Huynh, Q., Carruthers, T., Taylor, N.G. 2022. EcoTest: a proof of concept for evaluating ecological indicators in multispecies fisheries, with the Atlantic longline fishery case study. SCRS/2022/106.

Juan-Jordá, M.J., Murua, H., Arrizabalaga, H., Dulvy, N.K., and Restrepo, V. 2018. Report card on ecosystem-based fisheries management in tuna regional fisheries management organizations. Fish Fish. 19: 321-339.

14 Acknowledgements

EcoTest is funded by ICCAT as part of a commitment to the Global Environmental Facility (GEF) Common Oceans ABNJ Tuna Project Fund.

Many thanks to The Ocean Foundation for funding Phases I and II of the EcoTest framework.

Many thanks to the ICCAT Ecosystem working group for comments on earlier versions of EcoTest: Nathan Taylor, Alex Hanke, Sachiko Tsuji, Guillermo Diaz, Diego Alvarez, Laurie Kell, Maria Juan Jorda, Andres Domingo.

Blue Matter: Tom Carruthers, Adrian Hordyk, Quang Huynh

OpenMSE was developed with support from the Natural Resources Defense Council (NRDC), the Gordon and Betty Moore Foundation, the Packard Foundation, the Marine Stewardship Council, Fisheries and Oceans Canada (DFO), the U.S. National Oceanic and Atmospheric Administration, the International Commission for the Conservation of Atlantic Tunas (ICCAT) and The Ocean Foundation.

OpenMSE continues to be developed with the support of the The Ocean Foundation.